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Abstract 


Employing  the  Galerkin  method,  we  find  altogether  four  solutions  for  the 
Navier-Stokes  equation  describing  the  airflow  around  a  fluid  sphere.  Two  solutions  are 
real,  and  two  are  complex.  Of  the  two  real  solutions,  one  is  a  standard  solution  described 
by  Kawaguti  some  time  ago.  A  new  real  solution  is  distinctly  different  from  the  standard 
one  and,  as  such,  gives  a  qualitatively  different  description  of  the  flow  around  a  sphere. 
For  large  Re5molds  numbers,  this  new  solution  should  be  appropriate  for  deducing  the 
critical  forces  on  the  fluid  sphere  responsible  for  its  breakup. 
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Executive  Summary 


In  support  of  the  Air  Defense  Intercept  effort,  the  author  proposes  to  model  the  fluid  breakup 
through  the  Navier-Stokes  equation  by  studying  the  flow  instability  aroxmd  the  spherical  fluid  at 
very  high  air  velocities.  Essentially,  the  air  treated  here  as  a  fluid  exhibits  highly  nonlinear 
behavior.  The  breakup  of  the  spherical  fluid  is  to  be  connected  to  secondary  flow  and  flow 
instability  of  the  air  around  the  spherical  fluid.  In  addition  to  the  standard  real  Galerkin  method 
solution  of  the  Navier-Stokes  equation,  valid  mostly  for  small  Reynolds  numbers,  a  new  real 
Galerkin  method  solution  has  also  been  found.  It  is  distinctly  different  from  the  standard  one 
and,  therefore,  should  be  capable  of  describing  the  flow  at  large  Reynolds  numbers.  As  such,  it 
is  appropriate  for  deducing  the  breakup  of  the  ejected  fluid. 
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1.  Introduction 


1.1  Purpose.  The  purpose  of  this  report  is  to  describe  and  illustrate  the  numerical  Galerkin 
method  solution  of  the  Navier-Stokes  equation  as  applied  to  the  breakup  of  an  agent  released 
from  a  missile. 

1.2  Background.  After  a  dense  fluid  is  ejected  by  the  impact  from  the  intercept,  it  travels  at 
very  high  (most  often,  sonic)  speeds  through  the  air.  In  what  follows,  the  author  assumes,  to  a 
good  approximation,  that  the  ejected  fluid  is  of  a  spherical  shape.  However,  the  secondary  flow 
and  flow  instabihty  of  the  air  (treated  here  as  a  fluid)  around  the  spherical  fluid  believed  to  be 
responsible  for  the  fluid  breakup  represents  a  highly  nonlinear  problem  when  being  tackled  by 
the  Navier-Stokes  equation.  The  flow  (velocity)  field  can  be  computed  in  one  of  the  two  ways: 
by  numerical  methods  seeking  a  solution  for  velocities  only  at  a  set  of  selected  grid  points  in  the 
flow  region,  or  by  the  Galerkin  method,  in  which  the  solution  of  the  flow  field  is  approximated 
by  a  continuous  function  expanded  in  terms  of  polynomials  with  arbitrary  coefficients. 
Regardless  of  whether  the  nmnerical  or  the  Galerkin  method  is  used,  the  fact  that  the  dense  fluid 
is  assumed  to  be  of  spherical  shape  should  be  helpful  when  formulating  the  problem  through  the 
Navier-Stokes  equation  for  the  motion  of  air  as  a  fluid  around  it. 

1.2.1  General  For  the  description  of  the  air  streaming  aroxmd  the  spherically  shaped 
high-density  fluid  (Figme  1),  the  form  of  the  Navier-Stokes  equation  as  originally  given  by 
Hughes  and  Gaylord  [1]  is  particularly  useful  Then,  choosing  the  origin  of  the  spherical 
coordinates  at  the  center  of  the  sphere,  an  assumption  can  be  made  that  far  away  from  the  sphere 
the  flow  is  of  a  constant  speed  U  along  the  z-axis.  It  is  easily  seen  that  the  z-axis  represents  axis 
of  axial  symmetry;  the  flow  (velocity)  field  is  independent  of  the  azimuthal  angle  (coordinate) 

In  essence,  there  are  only  two  independent  velocity  components — ^Ur,  the  radial  component  and 
ue,  the  angular  component.  For  such  an  axisymmetric  configuration  for  which  the  velocity  flow 
field  is  independent  of  (|),  a  Stake’s  scalar  stream  function  \|/  can  be  introduced,  which,  in  turn, 
determines  u^  and  u®. 
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Figure  1,  Representation  of  the  Coordinate  System  With  the  Origin  in  the  Center  of  a 
Fluid  Sphere  Traveling  With  Very  High  (Sonic)  Speeds  Through  the  Air. 

The  boundary  conditions  must  be  taken  into  account  properly  and  are  actually  quite  simple. 
The  velocity  components  u^  and  u^  must  be  both  zero  on  the  sphere,  and  the  flow  far  away  from 
the  sphere  must  be  uniform.  These  boundary  conditions  are  now  easily  translated  to  the  Stokes 
scalar  stream  function  y. 


The  Navier-Stokes  equation  in  the  form  of  Hughes  and  Gaylord  [  1  ]  is  now  rewritten  for  the 
Stokes  scalar  stream  function  h/.  Although,  instead  of  two  vector  field  components,  u^  and  uq, 
there  is  now  a  differential  equation  for  just  one  scalar  function  \|/.  Even  so,  this  is  still  a  very 
complicated  differential  equation,  quadratic  both  in  differential  operators  and  the  function  v|/. 
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Specializing  now  to  the  Galerkin  method,  y  must  be  expanded  in  terms  of  a  complete  set  of 
basic  functions  Fj(r,|),  where  r  is  the  radial  distance  and  |  =  cos^.  In  order  to  satisfy  the 
boundary  conditions,  the  simplest  approach  is  to  make  Fj(r,|)  a  linear  combination  of  the 
Legendre  polynomials,  P,  (|) ,  where  the  coefficients  in  the  expansion  (r)  are  still  functions  of 
the  radial  distance  r. 

Determining  the  coeflScients  aj(r)  is  actually  a  very  difficult  task,  except  in  simplified 
special  situations.  When  determining  these  coefficients  on  the  sphere,  the  Navier-Stokes 
equation  must  be  satisfied  exactly  for  approximated  vj/;  away  from  the  sphere,  there  is  generally 
an  error.  The  demand  that  the  distribution  of  this  error  throughout  the  flow  field  vanish,  at  least 
in  principle,  supplies  a  set  of  algebraic  equations.  Finding  the  solution  of  these  equations  is  a 
major  task  and  it  can  be  done;  the  solutions,  once  found,  determine  the  expansion  coefficients. 
This,  in  turn,  gives  the  approximated  y,  which  then  becomes  a  known  function  of  r  and  0. 

Having  obtained  the  stream  function  vj/,  many  properties  of  the  flow  around  the  sphere  can  be 
deduced  (i.e.,  the  flow  pattern  since  it  is  now  easy  to  obtain  the  velocity  components  Ur  and  ue,  as 
well  as  their  derivatives).  However,  to  find  out  imder  which  conditions  (e.g.,  at  which  air 
speeds)  the  sphere  starts  to  break  up,  the  drag  of  the  sphere  must  be  computed.  This  can  be  done 
with  the  help  of  the  normal  and  shear  stresses  which  are  given  in  terms  of  the  derivatives  of  air 
velocity  components  at  the  surface  of  the  sphere.  In  any  case,  the  drag  of  the  sphere  consists  of 
two  parts:  Dp,  the  pressure  drag  (also  known  as  the  form  drag)  caused  by  the  normal  stresses  and 
Ds,  the  skin  friction  caused  by  shear  stress  on  the  surface  of  the  sphere.  To  evaluate  Dp,  the 
pressure  which  can  be  obtained  in  principle  fi'om  the  Navier-Stokes  equation  in  the  Hughes  and 
Gaylord  form  [1]  must  be  known. 

1.2.2  Threat.  The  methodology  presented  in  this  report  addresses  the  effects  of  agent-fluid 
break-up  threats  on  soldiers  and  a  military  weapon  S5^tem  from  both  ballistic  and  cruise  missiles. 
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1.3  Scope.  The  scope  of  this  rq)ort  is  to  utilize  a  physics  approach  to  solving  and 
understanding  problems  of  fluid  breakup  associated  with  ballistic  and  cruise  missiles. 


2.  Theoretical  Formulation  of  the  Method 


It  can  be  assumed  that  for  a  hurling  airborne  spherical  fluid,  the  air  itself,  as  long  as  there  are 
no  additional  forces  acting  on  the  system,  to  a  good  approximation,  is  incompressible;  that  is,  its 
density  is  practically  constant  in  time,  which  means: 


Dt 


V. 


Owing  to  the  continuity  equation. 


1  Dp 
p  Dt 


+  V-U  =  0. 


(1) 


This  means  that  the  velocity  of  the  air-fluid  element  is  divergenceless: 

Vu  =  0.  (2) 

The  origin  of  the  spherical  coordinate  system  is  situated  at  the  center  of  the  fluid 

sphere,  whose  radius  is  a  (Figure  1).  Far  away,  the  flow  of  the  air  is  of  constant  speed  and,  by 
defimtion,  in  the  z-direction,  u  =  zU .  Furthermore,  the  axisymmetric  configuration  (the  velocity 
flow  field  of  the  air  around  the  fluid  sphere  is  independent  of  <}>)  is  exploited.  For  such  a  flow,  a 
Stokes  stream  scalar  fimction  ^  is  introduced,  which,  in  turn,  defines  the  r  and  ^components  of 
the  velocity  as 


4 


u 


(3) 


1  d\p  _  \  dyp 

sin^  dO  cosd 


and 


1 

r  sin  d  dt 


(4) 


In  what  follows,  the  form  of  the  incompressible  Navier-Stokes  equation  as  given  by  Hughes 
and  Gaylord  [1]  is  utilized  as  follows: 


— -ux(Vxu) 

at  ^  ’ 


(  \  \ 

=  -V  p+— pu^  -/iVx(Vxu), 

V  2  j 


(5) 


where  |i  is  the  kinematic  viscosity  of  the  air. 

As  mentioned  in  section  2,  a  very  important  and  rather  helpful  quantity  when  trying  to  solve 
equation  (5)  is  the  dimensionless  Reynolds  number  Re,  defined  generally  as 

Re  =  — » V-=vp,  (6) 

V 

where  U  and  L  denote  a  typical  flow  speed  and  a  characteristic  length  scale,  respectively  [2]. 
For  example,  for  air  flowing  around  the  dense  fluid  sphere,  L  =  2a,  where  a  is  the  radius  of  the 
sphere.  U,  of  course,  is  simply  the  magnitude  of  the  incoming  velocity  of  the  air  toward  the 
sphere.  Equation  (6)  shows  the  connection  between  the  ordinary  iji)  and  the  kinematic  ( p ) 
viscosity,  where  p  is  the  density  of  the  fluid,  the  air  in  our  case. 

Already  from  equation  (5),  the  usefulness  of  the  Reynolds  number  can  be  seen  by 
estimating  the  orders  of  magnitudes  for  inertia  and  viscous  terms.  Namely,  firom 
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-pux(Vxu)  =  -(p/2)Vu^  +p(u-V)u,  the  inertia  term  is  identified  with  p(u-V)u  while,  on 

the  other  hand,  from  —  pV  x  (V  x  u)  =  — /iV(V  •  u)  +  pV^u  =  pV^u,  the  viscous  term  is  identified 

with  pV^u .  Assuming  now  that  the  derivatives  of  the  velocity  components  typically  change  by 

amounts  of  order  U/L  over  distances  of  order  L ,  the  following  order  of  magnitude  estimates  for 
the  inertia  and  viscous  terms  in  equation  (5)  are  obtained: 

inertia  term:  p  |(u  •  V)u  |  =  p  0(U^  /L) ,  (7) 

and 

viscous  term:  |pV^u|  =  pO(U/L^) .  (8) 


Hence,  it  can  be  deduced  that 


inertia  term 
viscous  term! 


=  0 


v\HV 


=  0(RJ. 


(9) 


So,  it  can  be  seen  immediately  that  » 1  corresponds  to  the  motion  of  fluid  of  small 

viscosity  and/or  large  L.  However,  R  « 1  does  not  always  mean  the  motion  of  fluid  of  large 
viscosity;  R  « 1  can  also  be  achieved  with  small  L. 


With  the  assumption  of  steady  flow  (du/at  =  0 ),  after  taking  the  curl  of  both  sides  of  the 
equation  (variations  of  p  and  p  with  respect  to  space  and  time  are  assumed  to  be  negligible),  and 
taking  into  accormt  equations  (3)  and  (4),  one  obtains  the  following  equation  for  Tf/ : 


1  (d^pd  dp  d  ^  ^dp  2dp^^.,.  2  . 

-JT  .  — T - — — +  2cot^  — - —  0  1^-= — D  p, 

r  smdydd  dr  dr  dd  dr  r  d9 )  ^  R  ^ 


(10) 
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where  the  operator  is  defined  as 


_  2  sind  d  f  1  d  ^ 

dr^  r^  ddysinddd) 

As  usual,  the  Reynolds  niimber  based  on  the  diameter  of  the  sphere  is 

„  2pUa 


(11) 


(12) 


Now,  it  is  clear  that  the  velocity  of  the  air-fluid  element  depends  on  r  and  0;  as  such,  it 
should  vanish  at  the  surface  of  the  sphere.  Therefore, 

H^e)  =  0,  (13) 


and 


^(a,0)  =  O.  (14) 

or 

Furthermore,  away  fi'om  the  sphere  r  ^  oo ,  it  is  obtained  firom  equation  (10)  that 

DV(r5^)  =  0  asr->oo,  (15) 


implying  that 


^  (r,0)  =  ^r^  sin^  5  as  r  ^  oo. 


(16) 
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which  is  the  statement  that  the  flow  far  away  from  the  body  is  uniform;  U.  is  the  initial  velocity 
of  the  air  in  the  z-direction. 

Equations  (13),  (14),  and  (16)  have  to  be  satisfied  throughout  the  calculations;  in  fact,  these 
equations  practically  require  the  use  of  Legendre  polynomials  when  expanding  \}/  in  terms  of  a 
complete  set  of  basic  functions. 

Without  actually  solving  the  problem  numerically  and/or  analytically,  the  behavior  of  the 
fluid  sphere  in  the  air  at  very  high  velocities  can  be  addressed.  Namely,  at  small  Reynolds 
numbers,  ,  the  drag  of  the  sphere  is  dominated  by  the  skin  fiiction;  while  for  R^  >  90,  the 

pressure  drag  becomes  more  important.  Since  R,  increases  as  the  velocity  around  the  sphere 

increases,  this  means  that  the  pressure  is  rather  important  in  the  breakup  of  the  fluid  at  higher  air 
velocities.  Now  one  can  proceed  phenomenologically  to  determine  how  the  breakup  of  the  fluid 
sphere  is  related  to  drag,  pressure  stress,  shear  stress,  or  the  pressure  itself;  however,  this  may 
take  some  time  to  figure  out. 

3.  Implementation  of  the  Galerkin  Method  to  Solve  the 
Navier-Stokes  Equation  for  Finite  Reynolds  Numbers 


As  mentioned  earlier,  the  calculations  are  easier  to  track  if  we  use  the  new  variable 


I  =  cos^. 


(17) 


With  this  variable,  the  nonlinear  differential  equation  from  (10)  now  becomes 


—  DV-  I  ^  ^  ^  ^ 

Re  r^a^ar 


2k  d\P 
l-k^  dr 


2^ 

Tdk 


dV  =  o, 


(18) 


where  now 
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(19) 


and 


D' 


5'  1-1'  a' 

ar''"  r' 


(20) 


Furthermore, 

D'^^'  =  0  asr->oo 


yields 


(21) 


\P(t,0  =  y  r'(l  -  ^')  as  r  ->  00 . 


(22) 


Now,  the  Galerkin  method  demands  that  the  solution  to  equation  (18)  or  equivalently  to 
equation  (10)  be  approximated  by  a  polynomial  formed  from  a  set  of  some  basic  functions. 
Taking  into  account  equations  (13-16)  and  (22),  it  is  not  difficult  to  see  that  these  basic  functions 
should  involve  the  Legendre  polynomials,  defined  here  by  Rodrigues’  formula  [3]: 


PJI)  = 


1  d" 
2"’m!d^'" 


(I'-r. 


(23) 


A  few  quotes  follow; 


•  Po(^)  =  l- 

•  P,(^)  =  ^ 

.  P,(^)=i(3^'-l). 
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•  P,(?)  =  i(55=-35). 

•  P.®  =  -i(355'-305'+3). 

o 


The  basic  functions  Fj(r,^)  are  chosen  in  the  following  form  [1,4]: 

F,(r,a  =  f.(r)[P,.,({)-Pi.,({)],  (24) 

where,  due  to  orthogonality  of  the  Legendre  polynomials,  the  orthogonality  of  F;  (r,4)  is 


The  solution  is  simply  approximated  by  the  following  linear  combination: 


U 


(25) 


(26) 


As  already  argued  by  Kawaguti  [4],  it  should  be  possible  now  to  choose  the  coefiScients 
f;  (r)  in  such  a  way  that  while  the  solution  in  equation  (26)  is  only  approximate,  all  the  boundary 
conditions  are  satisfied  exactly. 

For  example,  fj(r)  =  r^/3  reproduces  \f/  in  form  (equation  [22]);  however,  other  terms  are 
clearly  needed,  as  otherwise  conditions  in  equations  (13)  and  (14)  cannot  be  satisfied  [4],  The 
simplest  condition  to  satisfy  equation  (13)  is  to  change  f,(r)  =  r^  /3  into  (rV3)(l-a/r),  which 
also  satisfies  conditions  in  equations  (13)  and  (22).  Unfortunately,  even  for  n  =  2 ,  the  condition 
in  equation  (14)  is  not  satisified;  more  terms  as  a  power  series  in  r"'  are  needed  [4]. 
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With  Kawaguti  [4],  n  =  2  is  chosen  because 


p««)-p,e)=|(i-{’) 

and 

p,({)-p,(f)=Y(i-{‘)- 

We  write  for  ^ 


■^1  .  -^2  .  ^3  .  -^4  1 

> 

) 

Ip  P=  P» 

B,  B.  B.l 

Ip  /j 

/ 

(27) 


(28) 


(29) 


and 

0--.  (30) 

a 

Equation  (29)  automatically  satisfies  the  boimdary  condition  in  equation  (22)  of  viniform  flow 
at  infinity.  The  two  conditions  in  equations  (13)  and  (14)  at  the  surface  require 


(31) 


(32) 


(33) 
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and 


EiBi=0.  (34) 

i=l 

For  later  use,  the  following  identities  are  from  equations  (32-35): 

2(2  +  i)Ai=0,  (35) 

i=l 

and 

4 

2](3+i)Bi=0.  (36) 

i=l 

The  next  step  is  to  substitute  equations  (29)  and  (30)  into  equation  (18),  starting  with  just 
A,,  Aj,  B,,  andBj  while  assuming,  at  least  temporarily,  that  the  rest  of  coefficients  are  zero. 
It  is  easily  seen  that  this  is  not  sufficient.  Equations  (31—36)  are  solved  exactly  to  yield 
Aj  =  —2,  Aj  =3/2,  and  Bj  =  B2  =  0  and,  as  such,  give  a  solution  for  ^ .  However,  when  this  is 
substituted  into  equation  (18),  the  error  is  just  too  big  and  impossible  to  eliminate  since  the 
coefficients  have  already  been  determined.  Thus,  it  is  best  to  assume  that  all  A’s  and  B’s  from 
equations  (29)  and  (30)  are  different  from  zero.  The  author  may,  however,  later  try  to  see  what 
happens  if  some  of  them  (e.g.,  A4  and  B4)  are  put  to  zero.  Substituting  with  these  coefficients  for 
yp  into  equation  (18),  an  error  denoted  [1,  2,  4]  as  [NS]  is  obtained  on  the  left-hand  side  of  the 
Navier-Stokes  equation: 

[ns]  =  (1  -  )[g  0  (p)  +  ?  G 1  (p)  +  2  (p)  +  ^"G  3  (p)  +  (1  -  3^'  )(E  0  (p)  +  .  (p))] ,  (37) 

denoting 
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A  R 

a„  =—  and  b„  =^. 

"  p"  "  p” 


Functions  in  equation  (37)  are  expressed  as  follows: 


Go0>)  =  4  +90a,)-(2a,  +5a,  +9a42'’i  >  ^9) 

P  Re  Vi=l  A 


-Sbj  -211>4)+f  +35aj  +73a<) 

G.(p)  =  A  .  ^  )  ,  (40) 

+  Zb|  K-3b,-7bJ 


^+j;ai  (-6b, +15b,+42bJ+  j;(2+i)a,  (-2b, +3b, +7bJ 

4  I  2  i_i  )  J 


+  3(2a2  +  5a3  +9a4]j 


.  (41) 


G3(p)  =  -^|^5](3+i)b,J(2b3  -3b3  -7bJ, 


Eo(p)  =  4-  fp^-i;iaj(2b,-3b3-7bj-(8a3+25a3  +  54aj^bj  ,  (43) 

P  V  i=l  /  V  i=l  /_ 


^i(P)  =  y  I  SbiJ(6b)  -15b3  -42bJ-|^X;ib(J(2b,  -3b3  -7b4) 
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First,  the  Navier-Stokes  equation  is  satisfied  exactly  on  the  sphere  by  the  approximation  for 
,  equations  (29)  and  (30).  Remember  that  when  r  ->  a ,  p  ->  1 ,  so  that  a;  A-  and  bj  ->  Bj 
which,  with  the  help  of  equations  (3 1—36),  demands: 

9A2  +  35A3  +  9OA4  =  0,  (45) 


and 


B, -6B3-2IB4  =0.  (46) 

With  Kawaguti  [4],  in  order  to  minimize  the  error,  [NS]  is  expanded  in  terms  Fi(p,^)  from 
equation  (24)  with  ff(p)  =  p  and  each  term  must  be  zero.  This  is  achieved  if  each  coefficient 
in  the  expansion  is  zero;  for  n  =  2 ,  these  conditions  are  obtained: 

00  1 

Jdpp-‘  Jd^(Po(0- P2(^))[NS](p,^)  ^  NSl  =  0,  (47) 

1  -1 


and 


CO  1 

Jdpp-^'  /d§(P.(4)-P3(4))[NS](p,^)  =  NS2  =  0.  (48) 

1  -1 

The  integrals  in  equations  (47)  and  (48),  although  tedious,  are  nevertheless  straightforward  to 
cany  out.  To  facilitate  this,  the  following  functions  are  defined: 

fi  (^)  =  (1  -  ^' ),  i  =  0,1,2,3;  (49) 


and 
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lii(e)  =  S'(l-i")(l-35'),  i=04. 

(50) 

Then  from  equations  (47)  and  (48),  it  is  easily  seen  that  in  variable  ^ ,  the  following  integrals  are 

needed: 

If,»  =  id?f,(l)(P,(9-P,(f)), 

-1 

(51) 

If,„  =  j'd|f,({)(P,«)-P,(|)),i  =  0,l,2^, 

-1 

(52) 

ni,»  =  jfdShiKKp.o-PjK)), 

-1 

(53) 

and 

ih,„ = idfh,({)(p, «)-?,({)), 1=0,1. 

-1 

(54) 

All  are  zero  except  for  the  following: 

Ifoo2  =8/5,  If3a2=  8/35, 

(55) 

Ifj,3  =8/21,If3,3  =8/63; 

(56) 

and 

Ihoo2  —  32/35 . 

(57) 
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Integrals  in  p  must  be  calculated.  They  are  defined  as  follows: 


W 

IGi,  =  Jdp  G,(p)p-\  i  =  0,l,2,3. 


(58) 


TO 

IG,,  =  jdpG,(p)p-\  i  =  0,1,2, 3; 


(59) 


lEj,  =  Jdp  Ej(p)p-‘,  i  =  0,l. 


(60) 


and 


IEi2  =  jdp  Ei(p)p■^  i  =  0,1 , 


(61) 


In  these  notations,  equations  for  NSlandNS2  (equations  [47]  and  [48],  respectively)  are 
now  written  as  follows: 


NSl  Ifoo2l^01  1^102^^11  ^202^^21  ^302^^31  ^002  ^01  ^102^11  ~  (^2) 

and 

NS2  =  Ifo,3lGo2  +  If„3lG,2  +  If2,3lG22  +  ^3,31032  +  Dio, 3^02  +  Dl„3lE,2  =  0  .  (63) 

These  integrals  from  equations  (52-61)  have  been  evaluated  using  the  third  edition  of 
Mathematica  [5],  yielding  finally  for  NSl  and  NS2  the  following: 
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NS1  = 


8[2880A,  +  7(-  45  +  1035A;  +  I96OA3  +  3204A4)]Bi  32X36^ 


11025 


+ 


46O8A4B2  8B 
385 
56A3B4 


35 


3  ,  64A,B3  64A3B3  232A4B3  ,  32B4  ,  176A,B4 

_i - j - 1-  - 
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35 


35 


35 


25 


I28A4B4  64A3  I44A4 

4  + - 2-  +  . - ^ 


(64) 


15  455  R,  R, 

16A2(-1386  +  176B3R,  -99B3R,  -444B4RJ  ^ 
1155R.  "  ’ 


64A^  4OA3'  I46A4  2336A,A4  1168A^  I6A3  ,  ^  x 

NS2  = - -+ - -  + - -  + - i— ^  + - -  + - ^(65  +  91A,  +2I6A4) 

35  9  21  231  147  273  ^  ' 


8A 


693 

21 


7  A p/-At.  /\'>p  A  \  8B?  128B]B2  16B,B3  64B2B3 

^(l32  +  176A,  +564A3  +935A4) - i - ^ - ^ 

^  \  1  3  4/  me  ^01 


32^B2  ^  I6B2B4  ^  352B3B4  ^  64B:  32Bi  ^  48B3  ^ 


63  567  105  231 

2  ooT>  /lotD  ^43 


(65) 


693 


27 


273 


63  21R. 


R. 


3R. 


=  0. 


Equations  (31-34),  (45),  (46),  (64),  and  (65)  are  the  eight  equations  necessary  to  solve  for  the 
coefficients  Aj,A2, A3,A4,B,,B2,B3,andB4 .  Unlike  equations  (31-36),  (45),  and  (46), 
which  are  linear,  equations  (64)  and  (65)  are  quadratic  in  these  coefficients;  this  fact  makes 
solving  these  equations  simultaneously  rather  difficult.  Nevertheless,  equations  (31-36),  (45), 
and  (46)  are  solved  expressing  six  of  the  coefficients  first  in  terms  of  Aj  and  B,  and  then,  as  an 
exercise,  a  different  set  of  six  coefficients  in  terms  of  A4  and  B4 .  The  results  are  as  follows: 


15 


1 


A,=  -—(S  +  5A,),  A3  =— (17  +  7A,),  A4  =-— (95  +  34A3), 


29 


29 


23  19  5 

B2-— B„B3-^B„B4--B,; 


(66) 


or 


A,  =-^(95+58A,),  A,  =^(7+10A,).  A,  =-^(3+14A,), 
34  34  34 


(67) 
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Utilizing  equation  (66)  and  substituting  it  into  expressions  of  NSl  and  NS2 ,  equations  (64)  and 
(65),  explicit  expressions  are  obtained  for  them  as  follows: 


NS1(A,,B,,R  401625 h-87302A,)B,  ^  72(9  +  2A,) 

'  137162025  29R, 


(68) 


and 


NS2(A.,B,,RJ= - i - 

64438719345R, 

[70916685840B,  +  3645(6413185 +  756034A,  +32880Af)R,  -24207344BfR  J, 


(69) 


which,  when  solved  formally,  yield 


NS1  =  0-  A  4725(162162 +  85B,  Re) 

'  - 1702701 00 +  87302B,  Re’ 


(70) 


NS2  =  0:  A,(±)  = - - 

239695200R, 

,  (71) 

-2755743930Rg  ± 7830^2^^ (- 26405533440B,  -2805569379Rj  +9013504BfR J2 

Unfortunately,  this  is  the  maximum  progress  achieved  without  specifying  the  values  for  the 
Reynolds  number  R^ .  Namely,  in  order  to  get  the  allowed  Bi’s,  the  Aj’s  from  equation  (69) 

must  be  equated  with  the  Ai  from  equation  (68).  The  resulting  equation  is  quadratic  for  each  of 
the  two  B]  s,  yielding  altogether  four  B] ’s  as  well  as  four  Ai ’s.  As  a  consequence,  there  are  four 
sets  of  coefficients:  A, ,  A^ ,  A3 ,  A^ ;  B, ,  B^ ,  B3 ,  B, .  However,  for  R,  between  0  and  1,000, 
only  two  sets  are  real,  while  the  other  two  sets  are  complex;  only  real  coefficients  can  define  the 
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Stokes  stream  function.  Furthermore,  since  only  one  set  of  these  coefficients  is  calculated  in  the 
literature  [1,  4]  and  two  independent  solutions  for  the  Stokes  stream  function  are  offered,  this 
approach  shows  definite  progress  in  this  respect.  Also,  the  single  set  solutions  for  the 
coefficients  differ  numerically  between  references  [1]  and  [4]  as  well  as  with  the  author’s  precise 
calculations,  which  employed  the  third  edition  of  Mathematica  package  [5]. 

Specifying  the  Reynolds  number  R^,  Table  1  presents  the  calculations  of  pairs  of  real 
Ai,  Bi,  onl>^,  the  rest  of  the  corresponding  real  coefficients  A2,A3,A4;B2,B3,B4can  then  be 
obtained  from  equation  (66). 


Table  1.  Calculated  Coefficients  Ai  and  Bi  for  Two  Numerical  Solutions,  1  and  n,  of  the 
Navier-Stokes  Equation  by  Galerkin  Method.  Here,  the  Reynolds  Number  Is 
Varied  Between  1  and  10,000 


1 

_ n _ 1 

Re 

Ai 

Bi 

Ai 

Bi 

1 

-4.49912 

-0.188993 

22.6882 

2931.63 

10 

^.41227 

-1.89842 

20.057 

309.868 

^.41547 

-3.84605 

16.645 

171.199 

-3.7452 

-5.87988 

126.176 

hhesihh 

-3.21602 

-8.00948 

12.6375 

103.969 

-1.98757 

-10.2042 

11.459 

90.7646 

60 

-1.98757 

-12.3967 

10.5744 

82.0239 

70 

-1.38407 

-14.507 

9.88776 

75.8178 

80 

,  -0.830554 

-16.4721 

-9.34022 

71.1879 

90 

-0.338676 

-18.2581 

8.52323 

67.7492 

100 

0.0911558 

-19.8578 

8.52323 

67.7492 

150 

1.53085 

-25.5462 

7.33534 

56.2674 

2.30045 

-28.8338 

6.6967 

52.0863 

3.46462 

-34.1918 

5.68153 

45.9191 

4.03814 

-37.0206 

5.14878 

42.8967 

1,000 

4.15159 

-37.5961 

5.04031 

42.2982 

10,000 

4.55591 

-39.6918 

4.64482 

40.1622 
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In  terms  of  these  coefficients,  expressions  for  velocity  components  can  also  be  obtained  now 
by  carrying  the  differentiation  in  equations  (3)  and  (4).  The  results  are  as  follows: 

Ur  =^{^[^4  +p(A3  +A2P  +  A,p')]cos^  +  [b4  +p(b3  +  B2P  +  BjP^ )][1  +  3cos20]},  (72) 


and 


u«  = 


U 

^6 


sin^ 


+  (4A4  +  p(3A3  +  2A2P  +  AjP^ )+  (4B4  +  p(3B3  +  2B2P  +  BjP^  ))cos0)sin5 


(73) 


Utihzing  equation  (66),  these  equations  can  be  rewritten  in  terms  of  just  Aj  andBj.  The 
result  is  as  follows: 


\l^  =  a'Usin"  0 


'-95-34A, 

9(l7  +  7A,) 

15(8 +  5Ai)  A,T 

2 

< 

SSp'* 

29p' 

29p^  p  J 

-L 

5B,  19B, 

23B,  B,  1 

cos  6 

1 

w  1 

9p'‘  ^  9p' 

9p^  p 

(74) 


Ur  = 


1 

261p* 


Ux 


-29B,(-l  +  p)^(-5  +  9p)  T 

+  9(-95  +  306p-240p='  +2A,(-H- p)"(-17  +  29p))cos5  > 
+  87Bj  (- 1  +  pY  (-  5  +  9p)cos^  6 


(75) 


and 


20 


'  u 

sin^ 


P- 


26lp^ 


^^9(- 190  + 459/0 -240/o"  +  A,(-68  +  189p-150p'  +29p'))- 
+  29B,(-20  +  57p-46p'  +9p')cosP 


sin^  d 


(76) 


It  has  been  known  for  quite  some  time  that  the  Navier-Stokes  equation  has  a  steady  solution 
with  the  standard  coefficients  from  Table  1  (denoted  as  I)  when  the  Reynolds  number  is  below 
the  critical  point.  In  references  [1]  and  [4],  the  critical  Reynolds  number  is  deduced  to  be  about 
100.  Beyond  this  number,  the  steady  flow  becomes  superimposed  with  secondary  flows  and/or 
flows  with  vortices  that  become  dominant  as  the  Reynolds  number  approaches  1,000,  as 
indicated  by  Chow  [6].  On  the  other  hand,  Kawaguti  [4]  claims  that  the  solution  for  \|/  (with  the 

standard  coefficients)  describes  a  steady  flow  only  for  10  <  R,  <  80 .  At  R^  =  100.  and  beyond, 

the  flow  becomes  unstable  and  periodic  and  transforms  itself  into  flow  with  von  Karman’s 
vortex-street  or  with  distorted  loops  of  vorticity  arranged  with  some  measure  of  symmetry. 

Unfortunately,  as  the  Reynolds  number  increases  further  (to  1,000  and  beyond),  the  standard 
solution  for  i/'  becomes  more  difficult  to  connect  to  experiments,  hi  fact,  Kawaguti  [4]  has 
found  that  even  when  the  Reynolds  number  is  increased  indefinitely,  the  Galerkin  method  with 
the  standard  coefficients  (I)  would  predict  a  steady  solution  for  the  Navier-Stokes  equation, 
although  it  could  not  be  foimd  experimentally.  For  this  reason,  it  is  necessary  to  study  the 
solutions  to  the  Navier-Stokes  equation,  with  die  Galerkin  method  utilizing  the  second 
coefficients  (II).  It  is  possible  that,  since  coefficients  I  and  n  are  distinctly  different,  the  physics 
of  the  flow  will  also  be  different.  Hence,  more  than  likely,  the  solution  of  the  Navier-Stokes 
equation  with  the  coefficients  n  will  describe  the  turbulent  flow  which,  in  turn,  should 
correspond  to  the  fluid  breakup. 

Finally,  the  stability  of  the  solutions,  either  with  the  coefficients  I  or  n,  should  be 
investigated  by  the  perturbation  method.  Here,  following  Kawaguti  [4],  an  exponential  in  time 


21 


perturbation  is  added  to  the  steady  solution,  and  whether  this  perturbation  grows  or  disappears  is 
studied.  If  it  grows,  then  an  unstable  solution  indicates  the  likely  fluid  breakup. 

4.  Examples  of  Applications  of  the  New  and  Standard 
Solutions  of  the  Navier-Stokes  Equation 

Concerning  the  type  I  coefficients  entering  into  the  Galerkin  method  solution  for  the 
Navier-Stokes  equation,  by  using  a  special  type  of  perturbation,  Kawaguti  [4]  deduced  that  the 
critical  Reynolds  number  (below  which  the  airflow  around  the  spherical  fluid  is  steady)  is  about 
51.  Experimentally,  however,  it  is  known  to  be  around  100  and  beyond,  depending  on  the  fluid. 
Also,  the  type  I  coefficients,  while  reproducing  the  calculated  drag  in  a  fair  agreement  with 
experiments,  fail  to  do  so  for  vorticity  and  pressure  distributions  over  the  spherical  fluid  surface 
[4]. 


These  facts  alone  should  be  enough  to  pursue  the  solution  for  the  Navier-Stokes  equation 
using  the  Galerkin  method  with  the  type  n  coefficients. 

Figures  2-7  describe  examples  of  angular  dependencies  of  u^,  andu^  at  R,  =1,000, 
r  =  2a,  a  =  lm,  and  U  =  30ms  with  type  I  and  n  coefficients.  In  specific  evaluations,  when 
determining  R^,  the  viscosity  p  and  the  air  density  p  would  have  to  be  taken  explicitly  into 

account.  Nevertheless,  it  is  easily  seen  fi’om  these  examples  that  graphs  with  coefficients  n 
describe  a  different  physical  situation  than  the  graphs  with  the  coefficients  1.  Hence,  the  author 
believes  that  the  solution  with  the  type  n  coefficients  could  very  well  be  more  appropriate  to 
describe  the  fluid  breakup  at  very  high  (sonic)  velocities  of  the  fluid  in  the  air. 

5.  Conclusion  and  Recommendations 

Although  much  has  been  accomplished  in  this  report  when  evaluating  the  stream  functions 
and  velocity  components  around  the  spherical  fluid,  more  could  and  should  be  done.  In 
particular,  evaluations  of  the  drag  coefficients  of  the  sphere,  as  well  as  the  vorticity  and  pressure 
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Figure  2.  Example  of  the  Angular  Dependence  of  the  Stokes  Stream  Function  \|/  Around 
the  Spherical  Fluid.  It  Was  Evaluated  With  Coefficients  I  at  =  1,000,  r  =  2a , 
a  =  1  m,  and  U  =  30  ms"' . 


ue 


Figure  4.  Example  of  the  Angular  Dependence  of  the  Angular  Component  of  the  Air 
Element  Velocity  Ug  Around  the  Spherical  Fluid.  It  Was  Evaluated  With 
Coefficients  I  at  =  1,000,  r  =  2a ,  a  =  1  m,  and  U  =  30ms“' . 

w 


Figure  5.  Example  of  the  Angular  Dependence  of  the  Stokes  Stream  Function  v}/  Around 
the  Spherical  Fluid.  It  Was  Evaluated  With  Coefficients  II  at  R,  =  1,000 , 
r  =  2a ,  a  =  1  m,  and  U  =  30  ms  ' . 
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Figure  6.  Example  of  the  Angular  Dependence  of  the  Radial  Component  of  the  Air 
Element  Velocity  u^  Around  the  Spherical  Fluid.  It  Was  Evaluated  With 
Coefficients  II  at  R  =1,000,  r  =  2a,  a  =  1  m,  and  U  =  30ms"' . 


distributions  over  the  spherical  surface,,  should  be  done.  These  would  be  needed  not  only  to 
make  intelligible  comparisons  with  experiments  but  also  to  develop  a  full  picture  of  the  fluid 
breakup  after  leaving  either  a  ballistic  or  cruise  missile. 
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